Chapter 13 Source of functional differences

# A tibble: 3 × 2
  type_river count
  <chr>      <int>
1 core          25
2 endemic        7
3 marginal      11
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************

13.1 Taxonomic variation

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  group_by(type,sample) %>% 
  summarise(fraction=sum(count)) %>% 
  group_by(type) %>% 
  summarise(mean(fraction))
`summarise()` has grouped output by 'type'. You can override using the `.groups` argument.
# A tibble: 3 × 2
  type     `mean(fraction)`
  <chr>               <dbl>
1 core               0.901 
2 endemic            0.132 
3 marginal           0.0874

13.2 Genome size

# A tibble: 3 × 2
  type     genome_size
  <chr>          <dbl>
1 core        3337341.
2 endemic     3753149.
3 marginal    3281875.
# A tibble: 3 × 9
  .y.    group1  group2      n1    n2 statistic     p p.adj p.adj.signif
* <chr>  <chr>   <chr>    <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 length core    endemic     25     7        69 0.42  0.973 ns          
2 length core    marginal    25    11       139 0.973 0.973 ns          
3 length endemic marginal     7    11        44 0.659 0.973 ns          

## Functional differences between fractions

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>%
    ggplot(aes(x=type, y=value, group=type))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type, p.adjust.method = "BH")
# A tibble: 3 × 10
  trait .y.   group1  group2      n1    n2 statistic     p p.adj p.adj.signif
* <chr> <chr> <chr>   <chr>    <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 mci   value core    endemic     23     7        72 0.701 0.701 ns          
2 mci   value core    marginal    23    10       177 0.014 0.042 *           
3 mci   value endemic marginal     7    10        53 0.088 0.132 ns          
genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type),
           trait %in% c("B01","B02","B03","B04","B06","B07","B08","D01","D02","D03","D05","D06","D07","D09")) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type, p.adjust.method = "BH") %>% 
    print(n=100)
# A tibble: 42 × 10
   trait .y.   group1  group2      n1    n2 statistic     p p.adj p.adj.signif
 * <chr> <chr> <chr>   <chr>    <int> <int>     <dbl> <dbl> <dbl> <chr>       
 1 B01   value core    endemic     23     7      57.5 0.27  0.27  ns          
 2 B01   value core    marginal    23    10     168.  0.038 0.057 ns          
 3 B01   value endemic marginal     7    10      61   0.01  0.029 *           
 4 B02   value core    endemic     23     7      69   0.598 0.598 ns          
 5 B02   value core    marginal    23    10     173   0.024 0.073 ns          
 6 B02   value endemic marginal     7    10      53   0.088 0.131 ns          
 7 B03   value core    endemic     23     7      88   0.73  0.73  ns          
 8 B03   value core    marginal    23    10     160.  0.08  0.24  ns          
 9 B03   value endemic marginal     7    10      44   0.399 0.599 ns          
10 B04   value core    endemic     23     7      92.5 0.572 0.623 ns          
11 B04   value core    marginal    23    10     150.  0.182 0.546 ns          
12 B04   value endemic marginal     7    10      40.5 0.623 0.623 ns          
13 B06   value core    endemic     23     7      79   0.961 0.961 ns          
14 B06   value core    marginal    23    10     153   0.142 0.256 ns          
15 B06   value endemic marginal     7    10      49.5 0.171 0.256 ns          
16 B07   value core    endemic     23     7      65   0.471 0.471 ns          
17 B07   value core    marginal    23    10     174   0.022 0.042 *           
18 B07   value endemic marginal     7    10      58   0.028 0.042 *           
19 B08   value core    endemic     23     7      54   0.202 0.202 ns          
20 B08   value core    marginal    23    10     171   0.029 0.044 *           
21 B08   value endemic marginal     7    10      62   0.01  0.029 *           
22 D01   value core    endemic     23     7      98.5 0.39  0.585 ns          
23 D01   value core    marginal    23    10     172.  0.028 0.084 ns          
24 D01   value endemic marginal     7    10      39   0.73  0.73  ns          
25 D02   value core    endemic     23     7      66.5 0.508 0.508 ns          
26 D02   value core    marginal    23    10     168.  0.038 0.113 ns          
27 D02   value endemic marginal     7    10      53   0.086 0.13  ns          
28 D03   value core    endemic     23     7      84   0.883 0.883 ns          
29 D03   value core    marginal    23    10     182.  0.01  0.029 *           
30 D03   value endemic marginal     7    10      56.5 0.038 0.057 ns          
31 D05   value core    endemic     23     7      90   0.666 0.666 ns          
32 D05   value core    marginal    23    10     169   0.034 0.104 ns          
33 D05   value endemic marginal     7    10      46   0.315 0.473 ns          
34 D06   value core    endemic     23     7      64.5 0.447 0.447 ns          
35 D06   value core    marginal    23    10     163   0.062 0.092 ns          
36 D06   value endemic marginal     7    10      57   0.034 0.092 ns          
37 D07   value core    endemic     23     7      80.5 1     1     ns          
38 D07   value core    marginal    23    10     180.  0.012 0.036 *           
39 D07   value endemic marginal     7    10      51.5 0.11  0.165 ns          
40 D09   value core    endemic     23     7      93   0.537 0.806 ns          
41 D09   value core    marginal    23    10     136.  0.387 0.806 ns          
42 D09   value endemic marginal     7    10      35   1     1     ns          
genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type),by="genome") %>%
    pivot_longer(!c(genome,type),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type)) %>% 
    ggplot(aes(x=type, y=value, group=type))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

13.3 Functional differences between high and low endemisms

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment)) %>%
    ggplot(aes(x=type_prevalence_environment, y=value, group=type_prevalence_environment))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    to.domains(., GIFT_db) %>%
    as.data.frame() %>%
    mutate(mci=(Biosynthesis+Degradation)/2) %>% 
    rownames_to_column(var="genome") %>% 
    select(genome,mci) %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment)) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type_prevalence_environment, p.adjust.method = "BH")
# A tibble: 1 × 10
  trait .y.   group1 group2    n1    n2 statistic     p p.adj p.adj.signif
* <chr> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 mci   value high   low        5     2         2 0.381 0.381 ns          
genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment),
           trait %in% c("B01","B02","B03","B04","B06","B07","B08","D01","D02","D03","D05","D06","D07","D09")) %>%
    ggplot(aes(x=type_prevalence_environment, y=value, group=type_prevalence_environment))+
      geom_boxplot()+
      facet_grid(~trait,  scales="free")

genome_gifts %>%
    to.elements(., GIFT_db) %>%
    to.functions(., GIFT_db) %>%
    as.data.frame() %>%
    rownames_to_column(var="genome") %>% 
    left_join(core_microbiota %>% select(genome,type_prevalence_environment, type),by="genome") %>%
    filter(type == "endemic") %>% 
    select(-type) %>% 
    pivot_longer(!c(genome,type_prevalence_environment),names_to = "trait",values_to = "value") %>% 
    filter(!is.na(type_prevalence_environment),
           trait %in% c("B01","B02","B03","B04","B06","B07","B08","D01","D02","D03","D05","D06","D07","D09")) %>%
    group_by(trait) %>%
    pairwise_wilcox_test(value ~ type_prevalence_environment, p.adjust.method = "BH") %>% 
    print(n=100)
# A tibble: 14 × 10
   trait .y.   group1 group2    n1    n2 statistic     p p.adj p.adj.signif
 * <chr> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>       
 1 B01   value high   low        5     2       2   0.381 0.381 ns          
 2 B02   value high   low        5     2       2   0.381 0.381 ns          
 3 B03   value high   low        5     2       5   1     1     ns          
 4 B04   value high   low        5     2       4   0.857 0.857 ns          
 5 B06   value high   low        5     2       1.5 0.241 0.241 ns          
 6 B07   value high   low        5     2       1   0.19  0.19  ns          
 7 B08   value high   low        5     2       2   0.329 0.329 ns          
 8 D01   value high   low        5     2       1   0.171 0.171 ns          
 9 D02   value high   low        5     2       0   0.079 0.079 ns          
10 D03   value high   low        5     2       0   0.079 0.079 ns          
11 D05   value high   low        5     2       1   0.19  0.19  ns          
12 D06   value high   low        5     2       2   0.381 0.381 ns          
13 D07   value high   low        5     2       2   0.381 0.381 ns          
14 D09   value high   low        5     2       2   0.285 0.285 ns          

13.4 MAGs in different locations and shared among locations

13.4.1 Core

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  filter(type=="core") %>% 
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    facet_nested(~factor(environment, labels=c("low" = "Low altitude", "high" = "High altitude")),  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12, lineheight = 0.6,face="bold"),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

13.4.2 Endemic

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  filter(type=="endemic") %>% 
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    ylim(0,1)+
    facet_nested(~factor(environment, labels=c("low" = "Low altitude", "high" = "High altitude")),  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12, lineheight = 0.6,face="bold"),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

13.4.3 Marginal

genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == sample)) %>% #append sample metadata
  filter(count > 0) %>% #filter 0 counts
  left_join(core_microbiota,by="genome") %>% 
  filter(type=="marginal") %>% 
  ggplot(., aes(x=sample,y=count, fill=phylum, group=phylum)) + #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1) + #plot stacked bars with white borders
    ylim(0,1)+
    scale_fill_manual(values=phylum_colors) +
    facet_nested(~factor(environment, labels=c("low" = "Low altitude", "high" = "High altitude")),  scales="free") + #facet per day and treatment
    guides(fill = guide_legend(ncol = 1)) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          strip.background = element_rect(fill = "white"),
          strip.text = element_text(size = 12, lineheight = 0.6,face="bold"),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black")) +
   labs(fill="Phylum",y = "Relative abundance",x="Samples")

13.4.4 Ordination

gift_pcoa <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame() %>%
    vegdist(method="euclidean") %>%
    pcoa()

gift_pcoa_rel_eigen <- gift_pcoa$values$Relative_eig[1:10]


# Get genome positions
gift_pcoa_vectors <- gift_pcoa$vectors %>% #extract vectors
  as.data.frame() %>% 
  select(Axis.1,Axis.2) # keep the first 2 axes

gift_pcoa_eigenvalues <- gift_pcoa$values$Eigenvalues[c(1,2)]

gift_pcoa_gifts <- cov(genome_gifts, scale(gift_pcoa_vectors)) %*% diag((gift_pcoa_eigenvalues/(nrow(genome_gifts)-1))^(-0.5)) %>%
  as.data.frame() %>% 
  rename(Axis.1=1,Axis.2=2) %>% 
  rownames_to_column(var="label") %>% 
  #get function summary vectors
  mutate(func=substr(label,1,3)) %>% 
  group_by(func) %>% 
  summarise(Axis.1=mean(Axis.1),
            Axis.2=mean(Axis.2)) %>% 
  rename(label=func) %>% 
  filter(!label %in% c("S01","S02","S03"))
set.seed(101)
scale <- 20 # scale for vector loadings
gift_pcoa_vectors %>% 
  rownames_to_column(var="genome") %>% 
  left_join(genome_metadata, by="genome") %>%
  ggplot() +
      #genome positions
      scale_color_manual(values=phylum_colors)+
      geom_point(aes(x=Axis.1,y=Axis.2, color=phylum, size=length), 
                 alpha=0.9, shape=16) +
      #scale_color_manual(values=phylum_colors) +
      scale_size_continuous(range = c(0.1,5)) +
      #loading positions
      geom_segment(data=gift_pcoa_gifts, 
                   aes(x=0, y=0, xend=Axis.1 * scale, yend=Axis.2 * scale),
                    arrow = arrow(length = unit(0.3, "cm"), 
                    type = "open", 
                    angle = 25),
                    linewidth = 0.5, 
                    color = "black") +
     #Primary and secondary scale adjustments
     scale_x_continuous(name = paste0("PCoA1 (",round(gift_pcoa_rel_eigen[1]*100, digits = 2), " %)"),
                      sec.axis = sec_axis(~ . / scale, name = "Loadings on PCoA1")
            ) +
     scale_y_continuous(name = paste0("PCoA2 (",round(gift_pcoa_rel_eigen[2]*100, digits = 2), " %)"),
                      sec.axis = sec_axis(~ . / scale, name = "Loadings on PCoA2")
            ) +
    geom_label_repel(data = gift_pcoa_gifts,
                     aes(label = label, x = Axis.1 * scale, y = Axis.2 * scale),
                     segment.color = 'transparent') +
    xlim(-4.5,4.5) + 
    ylim(-3,2.5) +
    theme_minimal() + 
    theme(legend.position = "none")

structural <- core_microbiota %>% 
  mutate(type = case_when(
      (type_prevalence == "endemic") & (type_prevalence_environment == "high") ~ "present_high",
       (type_prevalence == "endemic") & (type_prevalence_environment == "low") ~ "present_low")) %>% 
    filter(!is.na(type)) %>% 
    select(genome,type)

enriched <- ancom_mag_skin_table %>% 
  mutate(type = case_when(
      lfc_environmentlow < 0 ~ "enriched_low",
      lfc_environmentlow >= 0 ~ "enriched_high")) %>% 
  select(genome,type)


gift_pcoa_vectors %>% 
  rownames_to_column(var="genome") %>% 
  left_join(bind_rows(structural,enriched), by="genome") %>%
  left_join(genome_metadata, by="genome") %>%
  mutate(type=ifelse(is.na(type),"neutral",type)) %>%
  group_by(type) %>%
  mutate(x_cen = median(Axis.1, na.rm = TRUE)) %>%
  mutate(y_cen = median(Axis.2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot() +
      scale_size_continuous(range = c(0.1,5)) +
      geom_point(aes(x=Axis.1,y=Axis.2, color=type, size=length), alpha=0.9, shape=16) +
      geom_segment(aes(x = x_cen, y = y_cen, xend = Axis.1, yend = Axis.2, color=type), alpha = 0.5) +
      scale_color_manual(values=c("#5A99D2","#F16144","#B7B7B7","#225072","#8C2C1C"))+
    theme_minimal()

scale <- 15 # scale for vector loadings
gift_pcoa_vectors %>% 
  rownames_to_column(var="genome") %>% 
  left_join(core_microbiota, by="genome") %>%
  left_join(genome_metadata, by="genome") %>%
  filter(type != "absent") %>% 
  mutate(shape = case_when(
      type == "core" ~ "core",
      type == "marginal" ~ "marginal",
      type == "endemic" & type_prevalence_environment == "high" ~ "endemic_high",
      type == "endemic" & type_prevalence_environment == "low" ~ "endemic_low"
    )) %>% 
  group_by(shape) %>%
  mutate(x_cen = median(Axis.1, na.rm = TRUE)) %>%
  mutate(y_cen = median(Axis.2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot() +
      #genome positions
      geom_segment(aes(x = x_cen, y = y_cen, xend = Axis.1, yend = Axis.2, color=shape), alpha = 0.5, show.legend = FALSE) +
      scale_color_manual(values=c("#9E8F71","#5A99D2","#F16144","#B7B7B7"))+
      theme_minimal() 

   # theme(legend.position = "none")
# A tibble: 33 × 5
   gift   high   low    p_value p_adjust
   <chr> <dbl> <dbl>      <dbl>    <dbl>
 1 B0104 0.681 0.824 0.00578    0.0279  
 2 B0106 0.894 0.747 0.00228    0.0123  
 3 B0208 0.476 0.720 0.000794   0.00631 
 4 B0213 0.582 0.795 0.00000629 0.000788
 5 B0220 0.179 0.311 0.00770    0.0359  
 6 B0302 0.233 0.457 0.00193    0.0108  
 7 B0307 0.370 0.201 0.000193   0.00260 
 8 B0309 0.134 0.269 0.000956   0.00679 
 9 B0402 0.222 0.357 0.000658   0.00631 
10 B0403 0.264 0.398 0.0116     0.0474  
# ℹ 23 more rows

term df SumOfSqs R2 statistic p.value
environment 1 14.464879 0.1785608 6.438403 0.001
river 2 8.130211 0.1003629 1.809402 0.097
Residual 26 58.413069 0.7210764 NA NA
Total 29 81.008159 1.0000000 NA NA

gift high low p_value p_adjust
B0104 0.67421663 1.0000000 0.0018754356 0.005672049
B0105 0.73759193 1.0000000 0.0018754356 0.005672049
B0106 0.80621259 1.0000000 0.0017459182 0.005672049
B0204 0.48537488 0.8606981 0.0145424059 0.028623148
B0206 0.61242517 1.0000000 0.0017937971 0.005672049
B0207 0.43143388 0.7100000 0.0019749754 0.005695278
B0209 0.79386259 1.0000000 0.0018754356 0.005672049
B0210 0.55247656 0.8606981 0.0145424059 0.028623148
B0213 0.61084133 0.8944683 0.0145424059 0.028623148
B0214 0.22485035 0.6834048 0.0219467710 0.041233327
B0218 0.46423032 0.7889365 0.0256749426 0.046819013
B0220 0.07385765 0.7171750 0.0010154015 0.005672049
B0601 0.54689170 0.9155746 0.0079058708 0.017824145
B0602 0.63919238 0.9155746 0.0145424059 0.028623148
B0603 0.66787867 1.0000000 0.0018754356 0.005672049
B0604 0.43445679 1.0000000 0.0018754356 0.005672049
B0605 0.00000000 0.4221269 0.0045906627 0.011161611
B0702 0.66854993 1.0000000 0.0018754356 0.005672049
B0703 0.07420061 0.6700000 0.0004222482 0.005672049
B0704 0.48294837 1.0000000 0.0018754356 0.005672049
B0706 0.41331224 1.0000000 0.0018754356 0.005672049
B0707 0.76551369 1.0000000 0.0018754356 0.005672049
B0710 0.06335476 0.1400000 0.0007537527 0.005672049
B0711 0.14734772 0.2664340 0.0010154015 0.005672049
B0712 0.71438862 0.3700429 0.0020470591 0.005768985
B0803 0.73937330 1.0000000 0.0018754356 0.005672049
B0804 0.52266792 0.8606981 0.0145424059 0.028623148
B1028 0.03147905 0.1400000 0.0014496951 0.005672049
B1029 0.00000000 0.1444683 0.0045906627 0.011161611
D0102 0.33798384 1.0000000 0.0018754356 0.005672049
D0103 0.16414075 0.6700000 0.0063405256 0.014834437
D0104 0.31230836 0.8000000 0.0019585157 0.005695278
D0201 0.06134920 0.8100000 0.0007255498 0.005672049
D0202 0.16943443 0.8037702 0.0010154015 0.005672049
D0203 0.22479340 0.4890978 0.0008505648 0.005672049
D0205 0.03972843 0.3057359 0.0010154015 0.005672049
D0206 0.00000000 0.6115146 0.0002918961 0.005672049
D0207 0.11256387 0.6128250 0.0010154015 0.005672049
D0208 0.13134920 0.4595489 0.0009808556 0.005672049
D0209 0.14134920 0.6079485 0.0010066900 0.005672049
D0210 0.19001248 0.2500000 0.0018754356 0.005672049
D0211 0.00000000 0.2195918 0.0045906627 0.011161611
D0213 0.00000000 0.3093019 0.0002918961 0.005672049
D0302 0.00000000 0.3300000 0.0001880890 0.005672049
D0306 0.11242517 0.5000000 0.0013515402 0.005672049
D0308 0.05621259 0.4422127 0.0007395593 0.005672049
D0309 0.25000000 0.6055317 0.0007325316 0.005672049
D0505 0.21728768 0.5264768 0.0145424059 0.028623148
D0506 0.16892979 0.5000000 0.0018099404 0.005672049
D0507 0.19124896 0.6273362 0.0010154015 0.005672049
D0508 0.08094613 0.2586895 0.0219467710 0.041233327
D0509 0.48650715 0.6444683 0.0223620729 0.041386523
D0511 0.32436673 1.0000000 0.0018754356 0.005672049
D0512 0.22485035 0.6834048 0.0219467710 0.041233327
D0513 0.06335476 0.3709022 0.0010154015 0.005672049
D0601 0.36577691 0.5000000 0.0018754356 0.005672049
D0604 0.00000000 0.4221269 0.0045906627 0.011161611
D0610 0.06917013 0.0000000 0.0051195000 0.012208038
D0611 0.11242517 0.5000000 0.0013515402 0.005672049
D0613 0.87641709 1.0000000 0.0126508551 0.027521158
D0704 0.56941341 1.0000000 0.0018754356 0.005672049
D0708 0.11242517 0.3944683 0.0120900251 0.026770770
D0804 0.67738433 0.0000000 0.0006075647 0.005672049
D0807 0.69205678 0.1544254 0.0020946080 0.005771809
D0813 0.00000000 0.2889365 0.0045906627 0.011161611
D0816 0.09668565 0.4820515 0.0065015311 0.014929442
D0901 0.00000000 0.4221269 0.0045906627 0.011161611
D0907 0.22485035 1.0000000 0.0013515402 0.005672049

term df SumOfSqs R2 statistic p.value
environment 1 49.92847 0.4416516 12.05864 0.001
river 2 13.43533 0.1188447 1.62244 0.187
Residual 12 49.68566 0.4395037 NA NA
Total 15 113.04945 1.0000000 NA NA

# A tibble: 0 × 5
# ℹ 5 variables: gift <chr>, high <dbl>, low <dbl>, p_value <dbl>, p_adjust <dbl>

term df SumOfSqs R2 statistic p.value
environment 1 12.004006 0.2843473 2.144208 0.106
river 1 7.818642 0.1852056 1.396600 0.318
Residual 4 22.393358 0.5304471 NA NA
Total 6 42.216006 1.0000000 NA NA